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Abstract 

We discuss a definition of robust dominant eigenvector of a family of stochastic matrices. Our focus is on 
application to ranking problems, where the proposed approach can be seen as a robust alternative to the standard 
PageRank technique. The robust eigenvector computation is reduced to a convex optimization problem. We also 
propose a simple algorithm for robust eigenvector approximation which can be viewed as a regularized power 
method with a special stopping rule. 



1 Introduction 

There are numerous problems which can be reformulated as finding a dominant eigenvector of a stochas- 
tic matrix — computing invariant probabilities in Markov chains, various problems in econometrics, 
sociometry, bibliometrics, sport etc., see an excellent survey on history of such applications 3J. Recently, 
ranking in the Web (PageRank) attracted a lot of attention [H |S] . 

Let us briefly remind this approach and discuss some problems which arise in relation to it. There 
are n pages (nodes of a graph) connected with outgoing and incoming links. Then Xi — score of i-th 
page of the Web, i = 1, n is defined as follows: Xi — X^eL x jl n v where rij is the number of outgoing 
links of page j and Li is the set of pages linked to page i. Thus the score vector x satisfies the equation 

Px = x (1) 

with Pij = 1/nj if page j cites page i and Pjj = otherwise. I.e. x is an eigenvector corresponding to 
eigenvalue 1 of the (column) stochastic matrix P of size n x n. 

In the example below n = 7, outgoing links are denoted by arrows. The corresponding matrix P and 




Fig. 1: Test example 
its dominant (or principal) eigenvector x are 



*A. Juditsky is with LJK, Universite J. Fourier, BP 51, 38041 Grenoble Cedex 09, France juditskyOimag.fr 
tB. Polyak is with Institute for Control Sciences, 65 Profsoyuznaya str., 117997 Moscow, Russia boris9ipu.ru 



/ 








1 /o 

1/3 















1/2 






















1/2 


1 





1/2 

























1 
















1/3 


1/2 































1 


V 








1/3 








1 





x = (0,0,0,0,0,0.5,0.5) T . 

We see that the scores of pages 6 and 7 equal 0.5 while all other pages have zero scores. Of course it 
contradicts common sense and our intuition about "true scores" . There are several problems related to 
the above definition: 

1. We have assumed that rij ^ for all j. Yet, in applications, there are nodes with no outgoing 
links — dangling nodes (e.g. pdf files). One can easily tackle this problem by taking rij = n for 
all dangling nodes. However more complicated problem of "traps" (absorbing nodes 6 and 7 in the 
above example) remains open. 

2. In the case of disconnected subgraphs (typical for real-life problems) x is not unique. 

3. Small changes in P may result in dramatic perturbations of a; — the computed scores are non- 
robust. 

4. Power method with iterations Xk+i — Pxk do not converge for cyclic matrices; for acyclic ones it 
may converge slowly due to the small gap 1 — | A2 1 (here A2 is the second eigenvalue of P). 

PageRank medicine for these problems goes back to the pioneering paper [3]: allow the i-the page to 
receive a part of its score "for free" - replace the dominant eigenvector x of P with the (unique) solution 
to the equation 

x a = aPx a + (1 — a)e, et — 1/n, i = 1, n. 
In other words, the "true" adjacency matrix P is replaced with 

M = aP + (1 - a)E, = 1/n, i,j = 1, ...,n, (2) 

so that x a = Mx a is the unique dominant eigenvector of stochastic matrix M. Unlike a dominant 
eigenvector of P, x a is robust to variations in links, with fast convergence of power iterations (as a k ). As 
a justification the model of surfer in the Web is proposed in [4]: a surfer follows links with probability 
a and jumps to an arbitrary page with probability I — a. Then xf is the (asymptotically) average time 
share he spends on page i. 

Here we propose a different modification of the original score definition ([I]). Note that in most 
of applications nominal matrices are known approximately and they are subject to perturbations, so 
robustness issues become important. Our goal is to provide a robust counterpart to standard techniques 
in the spirit of Robust Optimization framework [7]. In the hindsight, it shares many common features 
with with Robust Least Squares [BJ. To the best of our knowledge this approach to ranking problems is 
new. 

2 Problem formulation 

We say that matrix P £ M. nXn is stochastic (column-stochastic) if Pij > Vi, j, J2i — 1 ^3- The 
set of all stochastic matrices is denoted S. The celebrated Perron-Frobenius theorem states that there 
exist a dominant eigenvector x £ £ (here S = {u £ M. n \ ^ i Ui — 1, Uj > 0} being the standard simplex 
of R n ) of the "nominal" matrix P: 

Px — x. 

We are looking for a robust extension of the dominant eigenvector. Let V C S stand for the set of 
'perturbed stochastic matrices F = P + £, where £ is a perturbation. The function max^g-p \\Fx — x\\, 
where || • || is some norm, can be seen as a measure of "goodness" of a vector x as common dominant 
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eigenvector of the family V . We say that the vector x is a robust solution of the eigenvector problem on 
V if 

x G Argmin < max \\Fx — x\ 

Let us consider some choices of the parameters of the above definition - uncertainty set V and the norm 
|| • || - for the PageRank problem. Recall that our motivation is to "immunize" the score x against small 
perturbations of the adjacency matrix P. We may consider, that the uncertainty of the j-th column 
of P, in terms of the £i-norm of elements, is bounded by Sj. In other words, if is the j-th column 
of the perturbation matrix, then ||[£]j||i < £j, j = 1, ...,n. From now on we assume that £j > 0, 
j = 1, ...,n. Further, we may fix the "total uncertainty budget" e > ||£||i, where with some notational 
abuse, ||£|ji = With a "natural choice" of the norm || • || = || • || l7 the definition |3j) now reads 

x^~' G Argmin I max ||(P + £)x — x||i 

e T Kb"=0, mUi^ej, U\\i<e} (4) 

(here e. L — i = 1, ...,n). Let rij be the number of outgoing links on page j. If we set Ej < then 
F = P + £ is a stochastic matrix. 

A different robust eigenvector may be obtained if, instead of fixing a bound for the £i-norm of the 
perturbation, we bound its Frobenius norm: \\£\\f < £■ When combining it with the choice of the norm 
|| • || = || • ||2 in ([3]), we come to the definition 

£(2) g Argmin ■! max || (P + £)x — x\\2 

e T [^=0, ||[^||i<£i, ||el|F<e}. (5) 

An interesting simple version of the problem ^ may be obtained when substituting the "column-wise 
constraints" on £ with the requirement that P + £ is stochastic: 

xS F ' G Argmin < max || (P + t;)x — x\\2 

F + £ is stochastic, ||£||f < s}- (6) 

Obviously there is a vast choice of uncertainty sets and norms, which will result in different definitions 
of robust eigenvector. Three definitions above may be viewed as a starting point for this study. 




3 Theoretical analysis 

Observe that the min-max problems involved in the definitions Q - ^ of the robust eigenvector are 
not convex-concave. Indeed, if H is the corresponding set of perturbation matrices, the function 

<t>B{x) — ma? || (P + X ~ x \\ 

cannot be computed efficiently. However, one can easily upper bound the optimal values of these problems 
by those of simple convex problems as follows. 



Proposition 1: Let 

and let || • || = || • ||i. Then 



~! = U G e T [£\j = 0, lll^lli < e h lldli < e}, 



<t>3.Ax) < <fi(x), ¥>i(x) = \\Px - x||i + egi{x), 
where the convex function g\ is defined as 

gi(x)= mm {||u||oo + V" — \v\j\. 

j 
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The proofs of the statements of this section are postponed till the appendix. 

Note that cj\ is & norm. Further, when E\ — ... — s n and s — — is integer, gi(x) is simply the sum of 
s largest in magnitude elements of x. It is worth to mention that computing g\(x) requires 0(n) (up to 
a logarithmic factor in n) arithmetic operations. 

We have an analogous bound for function (f>~.(x), involved in the definition ([5]): 

Proposition 2: Let 

E 2 = {£ e R nxn \ e T [^ = 0, < ej, U\\f < e}, 

and let || • || = || • || 2 . Then 

<f>E 2 (x) < ip 2 (x), p 2 (x) = \\Px - x\\ 2 + eg 2 (x), 
where the convex function g 2 satisfies 

g 2 {x)= mm {||u|| 2 + Y] ^>|j }■ 

x=u+v (. * — ' £ J 

j 

What was said about the numerical cost of computing gi(x) may be repeated here: computing g 2 (x) 
requires 0(n) (up to "logarithmic factors") arithmetic operations. 
We have a particularly simple and intuitive bound for (JsJ) : 

Proposition 3: Let 

= {P + £ is stochastic, ||£||_f < el- 
Then 

<h. F {x) < <Pf(x), <Pf(x) = \\Px - x\\ 2 + e\\x\\ 2 . 

Note that the upper bound ip F (x) for (j)s F is not very conservative. Indeed, the corresponding worst-case 
perturbation is 

^ _ e(Px - x)x T 
\\Px — x\\ \\x\\ ' 

and for F = P + £* we have J2i Fij — 1; however the condition Fij > (required for F to be stochastic) 
may be violated. 

The results above can be summarized as follows: 
let || • || (x) and || • ||( 2 ) be some norms (in the above examples we have || • ||m set to l\- or ^-norm, and 
II ' 11(2) = 9ii 52 or || • || 2 ). With certain abuse of terminology, we refer to the vector 

x e Argmin {<p(x) = \\Px - x|| (1) + e||a;||( 2 )} (7) 

as (computable) robust dominant eigenvector of the corresponding family of matrices. 
Let us discuss some properties of the robust eigenvector. 

1. x coincides with x if P is regular and e is small enough. Recall that stochastic matrix P is called 
regular, if its largest eigenvalue Ai = 1 is simple, while all other eigenvalues lay inside the unit 
circle 1 — |Aj| > p, > 0, i = 2, ...,n. For instance, matrices with all positive entries are regular. 
Thus dominant eigenvector of a regular matrix is robust with respect to small perturbations. This 
is the case of matrix M in (J2J). 

2. In the case ||x||( 2 ) = || ■ H2 vector x is unique due to the strict convexity of ||a;|| 2 on E. 

3. For large e vector x is close to e € S, = n^ 1 , i = 1, ...,n. 
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4 Solving optimization problem ([7]) 



For irregular matrices or for matrices with small gap fi > the robust dominant vector x can not be 
found in explicit form. To compute it one should minimize the function ip(x) over standard simplex 
E. In the above examples these are well structured convex optimization problems, and as such they 
can be efficiently solved using available optimization software. For instance, medium-size problems (say, 
when n < l.e3 — l.e4) may be solved using interior-point methods. For large-scale problems one can use 
methods of mirror-descent family [51IH1E] and solve a saddle-point reformulation of as well as various 
randomized algorithms (see [10l [11] and references therein) . Finally, for huge-dimensional problems Yu. 
Nesterov |12j proposed recently special randomized algorithms which can work with sparse matrices and 
n = l.e6 — l.e9. We do not discuss these methods here. Instead, we propose a "fast and dirty" method 
for minimization of ([7]), which happens to be a close relative of the PageRank technique. 

The method starts the minimization process from the point X\, which is the minimizer of ||a:||(2) over 
E. Then we apply the power method with averaging to minimize the first term in ip(x): 

xi + Px x + ... + P k ~ 1 x 1 

Xk = , 

k 

(note that power method without averaging may be non convergent, for instance for a cyclic P). We 
have 



\Px k - Zfc||(l) = 



P k x x - x\\(v\ 



< 



c 



k - jfe' 

that is \\Px k - x k \\ {1) = 0(l/fc). 

To fix the ideas, let us consider the case of Proposition |3j where the norms || • ||(x) and || • Hp) are the 
Euclidean norms. We have x\ — e, and when writing the method in recurrent form we get 

x k+1 = (l - ^) Px k + ^e. (8) 

When denoting M k = ( 1 — k^jj -P+Fpi E, — 1/n, the iterative method reads x k +\ = M k x k , x% = 

e. We apply the method until (p k = ip(xk) starts to arise. Thus we arrive to 
Algorithm 1 
begin: X\ = e 

fc-th iteration: x k +\ — M k x k , 



stop: <p k+ i > ip k ,x = x k 

and x is the approximate solution. Note that the numerical complexity of one iteration is dominated 
with that of multiplication Px (and, of course, we should not write M k in dense form for sparse P). 
If we compare Algorithm 1 with PageRank algorithm: 

Xk+i = Mxk, M = aP + (1 — a)E, a = const, 

we conclude that our algorithm is PageRank with varying a, equipped with special stopping rule. The 
proposed method is also closely related to the the regularized algorithm of [13] . However, unlike the 
method of [T3], Algorithm 1 the stopping rule is related to the uncertainty level. 

Choice of e is an important issue for the described approach. One may consider the following heuristics 
for Web ranking problems: we may assume that the number nj of outgoing links of the j-th page is 
known with accuracy ±1 for qn pages, for some q < 1. Keeping in mind that the average rij = m = 20 
for the whole Web, and taking q — 0.5 we get for uncertainty in P, measured in Frobenius norm, 
e ~ y/qn/m ~ 0.03-y/n = 30 for n = 10 6 . In small problems (like the examples in this paper) the choice 
e ~ 1 appears to be convenient. 



5 



5 Numerical results 



5.1 Simulated graphs 

To illustrate an application of the robust eigenvector, proposed above, let us go back to the example in 
the introduction. We present of Figure [5] the scores computed for this example using three algorithms: 
1) PageRank algorithm with a = 0.85, 2) exact minimization of ifp with e = 1, 3) Algorithm 1 with 

6 = 1. 



V eigenvector 
□ 'robust', s= 1 

♦ iterative, stopped at 3-th step 

• 'google', a= 0,85 
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Fig. 2: Results for test example 

We observe that the scores for all three methods are close and differ dramatically from the dominant 
eigenvector of the nominal matrix. The robust approach generates the scores which are acceptable from 
the "common sense" point of view. It is worth to mentioning that Algorithm 1 requires only 4 iterations 

(ip 4 > (p 3 ,x = x 3 ). 

We applied the proposed approach to the simulated examples proposed in [TSj. The corresponding 
graphs may be constructed in arbitrary dimensions (from modest to huge ones), with the dominant 
eigenvectors of matrices P and M which can be computed explicitly. Furthermore, the corresponding 
matrix-vector product Px can be computed very efficiently without storing the matrix P itself. 

Consider the graph (Model 1) of Figure [3] The nodes are shown in circles, the references are denoted 
by arrows. To get rid of (n, n) as a dangling node we assume the ability of equally possible jump from 
this node to any other. There are N = n 2 nodes in the model. We are interested in finding the vector 
x = (xij) of scores. 

Taking arbitrary value of in, say X\\ = 1, we get the system of equations equivalent to ([T]): 
xn = xu = l^i-i,! + 1, i=2,...,n, x n = l, 

•^ij 2 ljj 2 — 1 ^ ' ^ 2, . . . , 77- 1 , 
X n j Xj n X n j—i -\- 2"^n—l,j ~t~ 1) 3 ^3 1; 
^™ ^n— l,n ^n,n— 1 "T" F 

These equations can be solved explicitly (from in we find 221 etc). Then we can proceed to standard 
normalization and get the scores x*, = ^ x ' 3 — . 

Similarly the scores for PageRank matrix M can be found from equations (it is convenient to take 
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Fig. 3: Model 1 



y n = a, then y nn = No): 

yii=yu = a(lyi-i,i + 1 )i i = 2,...,n, 

Vij = a{\Vi-i,j + \Vij-i + 1), j,i = 2,...,n-l, 

Vnj = Vjn = a{Vn,j-l + \Vn-\,j + 1), j = 2, 71 - 1, 

y nn = a(y n _i,„ + y n , n -i + 1) 

and normalization i/*„ = ^V^' 3 — . 

Note that even very small perturbations may strongly corrupt the graph on Figure [3] Thus, taking 
large penalty parameter e, as it was suggested in the discussion at the end of Section ^] would lead to 
complete "equalization" of the scores. That is why we choose small penalty coefficient in this example 
(specifically, it was taken s = 0.01). 

Let us denote x a = y*. We compare now the score vector x* for nominal P with the score x a by 
PageRank and the results obtained by Algorithm 1. On Figure [4] we present the scores of diagonal 
nodes, i.e. Xi t i,i = 1, . . . ,n for n = 200 (that is N — 40000) and a = 0.85. Figure [5] shows the scores 
x n i, i = 1, •• • ,n of the last row. From the geometry of our graph (Figure [3]) it is natural to expect the 
"true" scores Xij to increase with i,j and the highest scores should be achieved in south-east corner of the 
square. The score vector, obtained as dominant eigenvector of P, possesses these properties. PageRank 
solution and robust solution provide both strongly "equalized" ranks; the robust solution appears to 
better fit the expected one. 

Let us modify slightly the graph of Model 1: in the new Model 2 we link the node (n, n) only to the 
node (1,1). For this model the system of equations for which defines the nominal scores is similar to 
the system for Model 1, and the dominant eigenvector of P and PageRank vector x a can be computed 
explicitly. On the other hand, it is obvious that the power method would not converge for this model. 
Indeed if a^i = L^y = 0, {ij} ^ {11} then after 2n — 1 iterations we return at the same state, and 
matrix P is cyclic. 

The results for Model 2 with the same parameters (n — 200, e = 0.01) are presented at Figures [6] and 
[7] In this example, same as for Model 1, the score vector obtained using PageRank is close to the robust 
one. And may be robust scores are slightly closer to our expectations of "common sense" ranks. 



7 




Fig. 4: Diagonal scores 



5.2 "Real-life" data 

We have tested the proposed techniques on the experimental datasets, used in experiments on link 
analysis ranking algorithms described in [2]. We present here experimental results for adjacency matri- 
ces from two refined datasets: Computational Complexity (CompCom) and Movies, downloaded from 
http://www.cs.toronto.edu/tsap/experiments7l CompCom data represent 789 pages with totalize 1449 
links (with average number of links per page is equal to 1.8365). Movies dataset has 4754 notes and 
17198 (average number of links per page is 3.6176). We compare high accuracy solution to the problem 
ffih (with the norms ||-||(i) and ||-||(2) are Euclidean norms) obtained using MOSEK optimization software 
PJ with the scores by Algorithm 1 and those obtained by the PageRank algorithm. The same value of 
penalization parameter e = 1 in Q is used in all experiments. On Figure [8] we present the 20 largest 
elements of the score vectors for the CompCom data. In this experiment the optimal value of the problem 
([7]) is 0.0587. The corresponding value for solution supplied by Algorithm 1 (after 3 iterations) is 0.0756. 
It should be noted that the high scores x% 45 and x% 54 by PageRank are pure artifacts - they correspond 
to the pages which receive numerous links but only reference each other. The techniques proposed in 
this paper seem better handle this sort of problem than the classical PageRank algorithm. 

Figure [9] shows the 20 largest elements of the score vectors for the Movies data. In this case the best 
scores computed by 3 competitors are quite close. The optimal value of the problem Q is 0.0288; the 
approximate solution supplied by Algorithm 1 (after 4 iterations) attains the value 0.0379. 

6 Directions for future research 

1. Structured perturbations. It is often natural to assume that the perturbation £ of nominal matrix 
P possesses some structure. For instance, we can consider £ = MAN, where < e, and M,N 
are rectangular matrices of appropriate dimensions, or £ = J^tf-Af, where ||t||2 < e, with fixed 
matrices A4 € R" xn . 

2. Implicit uncertainty model. Family V can be defined implicitly. Suppose that we can sample 
randomly from T 3 , so that random realizations Fi G V, i = 1, N are available. One may try to 
find x which is approximate dominant eigenvector for all of them. For instance, one can fix e > 
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Fig. 5: Last row scores 



and look for solution of a system of convex inequalities 

\\FiX - x\\ < e, i = l,...N, xeS. 
Such system can be solved iteratively using the mirror-descent algorithm. This approach is close 

to mug. 

3. Computational research. The development of numerical methods tailored for optimization problem 
([7| and their comparison is of great interest. 

4. Real-life examples. We plan to test extensively the proposed approach on the large-scale Web rank 
problems data available from different Internet databases. 
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.1 Proof of Propositions [T] and [2] 

Note that 

(f> Sl (x) = max||(£ + P-i>||i 

= max max w T (£ + P — T)x. 
feHi ueK n ,||«|| 00 <i 

Then z = £ T u satisfies 



Wi = 



E 



%2 u i 



<EE^-i^- 
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We conclude that 



4>~A X ) < fi( x ) '■— max u T (P — I)x + egi(x) 

uGK" Jl-ixjloo^l 



= \\(P-I)x\\ 1 +eg 1 (x), (9) 

where 



gi(x) = e max z x 

zeR™, \\z\\oo<e, \z\ 3 <e 3 

= max z T x. (10) 

*eR", l|z|U<l, \z\i<e 3 le 



Note that ( 10 ) is a strictly feasible conic optimization problem, when dualizing the constraints we come 
to 

gi(x) = min {HU + "fN}- 

Along with (|9| this implies Proposition [I] □ 

The proof of Proposition [2] is completely analogous. 

.2 Proof of Proposition [3] 

We start with the following lemma, which is a modification of Theorem 3.1 in [6]. 
Lemma 1: For a e R n , b e R m , £ E R nxm it holds 

max ||a + #|| 2 = ||a|| 2 + e||6|| 2 . 

llsl|F<e 

Proof: Indeed, ||a + £6|| 2 < ||a|| 2 + e||&||2, while for a ^ 0, £ = £* = [mmn equality holds and \\C\\f = e 
(for a = one can take £* = p^, ||c|| 2 = 1, c G K™ arbitrary). □ 
We can now bound 

max II (P + £)x — x\\ 2 = max \\(P + £)x — x\\ 

?6S F '" p+aes,m\F<£ 

< max \\Px — x + £a;|| 2 = \\Px — x\\ 2 + e||x|| 2 , 

||C||f<e 

where the last equality is due to Lemma [l] □ 
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Fig. 7: Last row scores for the second model. 
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Fig. 8: 20 highest scores for CompCom matrix 
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Fig. 9: 20 highest scores for Movies matrix 
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